% (* ::Package:: *)

% Medium Stacks
I = eye (3109,3109);
%%  NOx NO3
NOx_county_storage_M = zeros(3109, 3109);

run Med_Reset
NOx_Med_Stack = 28761.72.*1.35.*(NOx_Cal.*Med_Stack {1,1});
p=11;
for n = 1:S; 
    sprintf('Medium Stacks: %d/%d, Loop: %d of 25.', n, S, p)
    Emission_Plus = (Med_Stack {4,1}(:,2)'+I (n,:));
    NOx (:,3) =((Emission_Plus)*((NOx_Med_Stack)))';
   
    run Nitrate_Sulfate_Ammonium_Marginal_New   
    run PM_25_Health

        Damages = [All_Mort{Cause,1}];
        NOx_M(n,1) = ((sum(sum(Damages)))) - B_25_Primary_MD;
        
        NOx_county_storage_M(n,:) = all_mort_damages_county' - all_mort_dam_cnty_base';
        clear all_mort_damages_county
    clear Damages
end
clear NOx_Med_Stack

%% Primary PM 25
PM_county_storage_M = zeros(3109, 3109);

run Med_Reset
PM_Med_Stack = 28761.72.*(PM25_Cal.*Med_Stack {2,1});
p=12;
for n = 1:S; 
    sprintf('Medium Stacks: %d/%d, Loop: %d of 25.', n, S, p)
    Emission_Plus = (Med_Stack {4,1}(:,4)'+I (n,:));
    PM_25_Primary (:,3) =((Emission_Plus)*(PM_Med_Stack))';
    
    run Nitrate_Sulfate_Ammonium_Marginal_New  
    run PM_25_Health

        Damages = [All_Mort{Cause,1}];
        PM_25_M(n,1) = ((sum(sum(Damages)))) - B_25_Primary_MD;
        
        PM_county_storage_M(n,:) = all_mort_damages_county' - all_mort_dam_cnty_base';
        clear all_mort_damages_county
    clear Damages
end
clear PM_Med_Stack

%% SO2_SO4
SO2_county_storage_M = zeros(3109, 3109);

run Med_Reset
SO2_Med_Stack = 28761.72.*1.5.*(SO2_Cal.*Med_Stack {3,1});
p=13;
for n = 1:S; 
    sprintf('Medium Stacks: %d/%d, Loop: %d of 25.', n, S, p)
    Emission_Plus = (Med_Stack {4,1}(:,5)'+I (n,:));
    SO2 (:,3) = ((Emission_Plus)*(SO2_Med_Stack))';
  
    run Nitrate_Sulfate_Ammonium_Marginal_New    
    run PM_25_Health

        Damages = [All_Mort{Cause,1}];
        SO2_M(n,1) = ((sum(sum(Damages)))) - B_25_Primary_MD;
        SO2_county_storage_M(n,:) = all_mort_damages_county' - all_mort_dam_cnty_base';
        clear all_mort_damages_county
    clear Damages
end
clear SO2_Med_Stack
%% NH3 NH4
NH3_county_storage_M = zeros(3109, 3109);

run Med_Reset
NH3_Med_Stack = 28761.72.*1.06.* ((NH4_Cal.*Med_Stack {5,1}));
p=14;
for n = 1:S; 
    sprintf('Medium Stacks: %d/%d, Loop: %d of 25.', n, S, p)
    Emission_Plus = (Med_Stack {4,1}(:,1)'+I (n,:));
    NH3 (:,3) = ((Emission_Plus)*(NH3_Med_Stack))';    
     
    run Nitrate_Sulfate_Ammonium_Marginal_New    
    run PM_25_Health

        Damages = [All_Mort{Cause,1}];
        NH3_M(n,1) = ((sum(sum(Damages)))) - B_25_Primary_MD;
        NH3_county_storage_M(n,:) = all_mort_damages_county' - all_mort_dam_cnty_base';
        clear all_mort_damages_county
    clear Damages
end
clear NH3_Med_Stack

%% VOC
VOC_county_storage_M = zeros(3109, 3109);

run Med_Reset
VOC_Med_Stack = Med_Stack {2,1};
p=15;
for n = 1:S; 
    sprintf('Medium Stacks: %d/%d, Loop: %d of 25.', n, S, p)
    Emission_Plus = (Med_Stack {4,1}(:,6)'+I (n,:));
    A_VOC (:,3) = (28761.72.*(Emission_Plus)*(VOC_Cal.*VOC_Med_Stack))';    
     
    run Nitrate_Sulfate_Ammonium_Marginal_New
    run PM_25_Health

        Damages = [All_Mort{Cause,1}];
        VOC_M(n,1) = ((sum(sum(Damages)))) - B_25_Primary_MD;
        VOC_county_storage_M(n,:) = all_mort_damages_county' - all_mort_dam_cnty_base';
        clear all_mort_damages_county
clear Damages
end
